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Abstract 

Background: It has been shown previously that glucocorticoids exert a dual mechanism of action, entailing 
cytotoxic, mitogenic as well as cell proliferative and anti-apoptotic responses, in a dose-dependent manner on 
CCRF-CEM cells at 72 h. Early gene expression response implies a dose-dependent dual mechanism of action of 
prednisolone too, something reflected on cell state upon 72 h of treatment. 

Methods: In this work, a generic, computational microarray data analysis framework is proposed, in order to 
examine the hypothesis, whether CCRF-CEM cells exhibit an intrinsic or acquired mechanism of resistance and 
investigate the molecular imprint of this, upon prednisolone treatment. The experimental design enables the 
examination of both the dose (0 nM, 10 nM, 22 uM, 700 uM) effect of glucocorticoid exposure and the dynamics 
(early and late, namely 4 h, 72 h) of the molecular response of the cells at the transcriptomic layer. 

Results: In this work, we demonstrated that CCRF-CEM cells may attain a mixed mechanism of response to 
glucocorticoids, however, with a clear preference towards an intrinsic mechanism of resistance. Specifically, at 4 h, 
prednisolone appeared to down-regulate apoptotic genes. Also, low and high prednisolone concentrations up- 
regulates genes related to metabolism and signal-transduction in both time points, thus favoring cell proliferative 
actions. In addition, regulation of NF-KB-related genes implies an inherent mechanism of resistance through the 
established link of NF-kB inflammatory role and GC-induced resistance. The analysis framework applied here 
highlights prednisolone-activated regulatory mechanisms through identification of early responding sets of genes. 
On the other hand, study of the prolonged exposure to glucocorticoids (72 h exposure) highlights the effect of 
homeostatic feedback mechanisms of the treated cells. 

Conclusions: Overall, it appears that CCRF-CEM cells in this study exhibit a diversified, combined pattern of 
intrinsic and acquired resistance to prednisolone, with a tendency towards inherent resistant characteristics, 
through activation of different molecular courses of action. 
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Background 

Resistance to glucocorticoids (GC) is considered to be 
one of the most important factors in the prognosis of 
leukemia [1,2]. In a previous study, it has been shown 
that when a resistant T-cell leukemia cell line (CCRF- 
CEM) is treated with prednisolone, the drug exerts a 
dual (biphasic) effect on these cells [3]. At low doses, 
prednisolone has a mitogenic/anti-apoptotic effect, 
whereas at higher doses it manifests a cytotoxic/mito- 
genic effect. Also, it has been shown that the actual 
underlying effect of prednisolone, either mitogenic or 
cytotoxic, becomes apparent at 72 h of prednisolone 
exposure, providing evidence for activation of a cellular, 
homeostatic, feedback mechanism at the transcriptional 
or translational layer (protein synthesis) [3]. 

In addition, it remains elusive whether cells possess 
inherent mechanisms inducing GC tolerance on them, 
or their responce upon GC treatment is one of gradual 
adjustment, meaning that originally sensitive cells 
become resistant. Thus, as glucocorticoid receptor regu- 
lates directly or indirectly several thousands of genes, 
this partly refers to activation of genes related to anti- 
apoptosis and mitogenesis. In this sense, those mechan- 
isms may possibly, through intricate, regulatory actions 
and cross-talks, confer to the induction of resistance in 
leukemic cells. Apoptosis evasion, or proliferation stimu- 
lation are two alternative mechanisms through which 
cells exhibit resistance. In the present work we refer to 
acute lymphoblastic leukemia (ALL), though glucocorti- 
coid treatment belongs to the first-line of medications 
against lymphoid malignancies in general [4,5]. There is 
adequate evidence supporting a far more intricate 
mechanism of resistance to glucocorticoids than mere 
down-regulation of steroid receptors. 

In this sense, several, possible resistance mechanisms 
of leukemic cells to glucocorticoid administration have 
been proposed, like the presence of somatic mutations 
on the GR gene that may lead to aberrant regulation of 
the receptor through intracellular signaling. Besides, sev- 
eral polymorphisms, but not somatic mutations, have 
been found in normal and ALL populations, not linked 
to resistance or sensitivity induction though, either in 
vivo or in vitro[6,7]. Other GC resistance scenarios are 
emphasizing in defects in intracellular signaling path- 
ways that involve interactions of GR with other 
sequence-specific transcription factors, such as AP-1 
and Nuclear Factor kappa-light-chain-enhancer of acti- 
vated B cells (NF-kB) [8]. In a normal cell, ligand-acti- 
vated GR may potentially interfere with transcription 
factor c-Jun or p65 NF-kB and thereby repress genes 
promoting cell proliferation and cell survival [6,9,10]. 
GR-dependent inhibition of the transcription factor p65 
NF-kB, plays a significant role in the manifestation of 



apoptotic and anti-apoptotic effects of GR in leukemia 
cells and has been identified as a pivotal component of 
the mechanism of cancer cell resistance to chemother- 
apy [9], Previous studies of GC effects on leukemia cells 
identified c-myc and cyclin D3 as early GR-regulated 
targets, in GC-sensitive cells [11]. Further studies 
showed that introduction of a conditionally expressed 
cyclin-dependent kinase inhibitor pl6 (INK4A) gene, 
sensitized GC-resistant leukemia cells, through induc- 
tion of cell cycle arrest [12]. 

Thus, pl6 inactivation may change GR levels, affecting 
GR-mediated gene regulation and resulting in resistance 
to GCs. For this purpose, the parental CCRF-CEM cell 
line was chosen as the system of study for the effects of 
prednisolone treatment, a T-cell leukemia cell line char- 
acterized by a mutation (L753F) on one GR gene allele 
that impairs ligand binding [13]. It is known that both 
the DNA and ligand binding domains of the GR are 
required in order to repress NF-kB transactivation [14]. 
Interestingly, concerning the question whether this muta- 
tion would affect GC resistance, it has been reported pre- 
viously that both the GC-resistant, as well as the GC- 
sensitive CCRF-CEM subclones, express heterogeneous 
populations of the GR (GR wt /GRL753F) [15,16]. The 
CCRF-CEM cell line has been reported to be resistant to 
GCs, presumably due to the accumulation of more resis- 
tant variants after long periods of prolonged culture [17]. 
In addition, utilization of an in vitro system provides 
reproducibility, an expedient system to systematically 
examine the impact of intracellular signals and at the 
same time minimize the effect of undesired crosstalks 
introduced by other in vwo-participating systems. 

A detailed molecular explanation of the intricate 
mechanisms, underlying the resistance phenotype to GC- 
induced apoptosis, remains elusive. The present work 
proposes a rational computational framework in order to 
aid the elucidation of the question whether the system 
under study, has intrinsic or acquired mechanisms of 
resistance. Our presumption is that the system in study 
possesses an intrinsic mechanism of resistance to gluco- 
corticoids i.e. prednisolone. Using the proposed compu- 
tational analysis workflow, we have analyzed microarray 
data from two time points (4 and 72 h treatment) and 
three different concentrations (10 nM, 22 uM and 700 
uM). For the 4 h time point, we used a 1.2 k platform, 
comprising of cancer specific genes, which has been 
reported and analyzed previously [3,18]. In order to 
expand our view of prednisolone effects on the cell line, 
we used a 4.8 k platform. Genes included in the 1.2 k 
platform are also represented in the 4.8 k platform. Data 
analysis was performed in order to find groups of genes 
associated with characteristics related to anti-apoptosis 
and apoptosis, cell cycle arrest, drug resistance etc. 
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Methods 

Data collection 

The CCRF-CEM cell line was obtained from the Eur- 
opean Collection of Cell Cultures (ECACC). Concentra- 
tions of prednisolone (Pharmacia) were: 0 uM (control), 
10 nM, 22 uM, and 700 uM. In general, all three predni- 
solone concentrations correspond to in vivo dosages 
administrated intravenously to children at ages between 
1 month and 12 years old [3]. Specifically, the 10 nM 
and 700 uM prednisolone concentrations were chosen 
as indicative of manifestation of specific phenotypic 
effects, i.e. anti-apoptosis accompanied with mitogenic 
effect and cytotoxicity accompanied with resistance, as 
observed by flow cytometry [3]. Moreover, the high con- 
centration (700 uM) used, is similar to concentrations 
used in different studies in CCRF-CEM cells [19] as well 
as primary cell cultures derived from childhood ALL 
patients [20]. The 22 uM prednisolone concentration 
was chosen as an intermediate concentration between 
the aforementioned two. 

RNA was isolated with Trizol (Invitrogen Inc.) accord- 
ing to the manufacturer's instructions. At least 40 ug of 
RNA from each sample was used. cDNA microarray 
chips from two platforms (1.2 k and 4.8 k) were 
obtained from TAKARA (IntelliGene® Human Cancer 
CHIP Version 4.0 and IntelliGene® II CHIP, respec- 
tively). Hybridization was performed with the CyScribe 
Post-Labeling kit (RPN5660, Amersham) as described by 
the manufacturer [3]. The experimental setups consisted 
of the five following pairs: control vs. 10 nM predniso- 
lone at 4 h (designated as '1'), 10 nM vs. 700 uM pre- 
dnisolone at 4 h (designated as '2'), control vs. 700 uM 
prednisolone at 4 h (designated as '3'), 22 uM vs. 700 
uM prednisolone at 72 h (designated as '4'), and control 
vs. 700 uM prednisolone at 72 h (designated as '5'). In 
this study, triplicate hybridizations are utilized for the 4 
h experiments. Experimental pairs were co-hybridized 
on the same slide, each stained with a different fluoro- 
phore. Fluorophores used were Cy3 and Cy5. Slides 
were scanned with the ScanArray 4000XL microarray 
scanner. Images were generated with ScanArray micro- 
array acquisition software (GSI Lumonics, USA). Image 
analysis was performed with the ImaGene* 6.0 software 
(Biodiscovery Inc., USA). The raw datasets have been 
deposited in NCBI's Gene Expression Omnibus (GEO), 
and are accessible through GEO Series accession num- 
ber [GEO: GSE28154]. 

Data preprocessing 

A common for both platforms data preprocessing stage, 
namely the median intensity value in each channel, was 
applied to the raw data. Specifically, the well performing 
robust version of the robust loess-based background 



correction (rLsBC) approach, as proposed by [21], was 
applied. rLsBC assumes that the background noise 
affects the spot intensities in a multiplicative manner 
[22]. Instead of using the measurements of the local 
(feature-related) background for the correction, rLsBC 
utilizes the regression estimate of the logarithmic back- 
ground distribution B ' according to the logarithmic 
foreground intensity F ' for each channel (7?: Red and 
G: Green). Thus, rLsBC provides a robust estimation of 
the channel-specific background noise, utilized to back- 
ground-correct the logarithmic foreground intensities: 

F R,G = pR,G _ B R,G 

where F R,G is the logarithmic background-corrected 
foreground intensity, and B^' G the robust estimate of 
background noise, for each channel. The absolute back- 
ground-corrected foreground intensity ff' G for each 
channel is then calculated as: 

= 

In order to reduce the complexity of the data set, we 
followed the replicate averaging approach proposed by 
[23]. In this approach, instead of estimating a constant c 
and utilizing it to adjust each of the individual replicate 
measurements, equivalently the replicates were averaged 
by taking their geometric mean, that is: 

rR,G _ 3 /fR,C rR,G rR,G 

la ~ yln ' Jr 2 'Jr 3 * 

where f^' G is the (background-corrected) foreground 
intensity of the replicate r b i = 1, 2, 3, and f^' G is the 
averaged foreground intensity across all replicates (hen- 
ceforth referred simply as signal intensity), for each 
channel (Red and Green). 

Since outliers can significantly influence one or more of 
the subsequent processing steps, extreme outlier values, 
that is signal intensities deviating more than 3 interquar- 
tile distances from the first or the third quartile [24], 
were identified and excluded in an iterative process. 

The signal intensities of each dataset were further nor- 
malized in order to mitigate the effect of extraneous, non- 
biological variation in the measured gene expression levels. 
The robust version of the intensity-dependent scatter-plot 
smoother loess [25,26] with a quadratic polynomial model 
was applied to the M-A scatter-plot [27], where M and A 
are the log-ratio and log-mean, respectively: 

M = F R — F G 
F R +P G 
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where F R and F° are the corresponding logarithmic 
signal intensities for each channel. The smoothing para- 
meter of the loess procedure used equals to 10%, which 
was considered appropriate for the relatively small num- 
ber of probes attached in the microarrays. 

In order to identify and remove dubious features, a fil- 
tering approach was followed based on (i) the loop- 
design [28-31] used for the experimental setups at 4 h 
(experiments '1', '2', and '3'), and (ii) the philosophy of 
the replicate filtering approach proposed by [32]. In par- 
ticular, the fold changes of experiment T' (control vs. 10 
nM prednisolone) and experiment '2' (10 nM vs. 700 
uM prednisolone) should roughly equal the fold change 
in experiment '3' (control vs. 700 uM prednisolone), 
that is: 

fR fR fR 

Ji ii _ J3 

A G 7 2 G "/ 3 G 

where f R,G are the signal intensity of the experiment 

i, i = 1, 2, 3, for each channel, or equivalently, the fol- 
lowing quantity should ideally be equal to zero: 

/f±Jl\ 
LdF = log 2 J = 

= Mj + Mj — M3 

where M; is the log-ratio of the experiment i, i = 1, 2, 3. 
Thus, we sought to filter out features whose LdF deviates 
greatly from this expected value of zero. We calculated 
the mean and standard deviation (SD) of LdF, and elimi- 
nated features whose LdF is greater than 2 SD from the 
mean. The selected number of standard deviations of the 
mean ensures that on the one hand the features used for 
the subsequent analysis roughly comply with this loop- 
design rule with relatively high confidence, while on the 
other hand selection thresholds are not too stringent, 
and thus reject potentially interesting genes. 

Data integration 

Before proceeding to the cross-platform analysis, since 
two different microarray platforms (IntelliGene® Human 
Cancer CHIP Version 4.0 and IntelliGene® II CHIP) 
were utilized in the present study, we chose to perform 
data integration at a lower level [33], instead of conduct- 
ing a meta-analysis. That is, the preprocessed datasets of 
each array were combined to a single, unified dataset, in 
which standard statistical procedures were applied. 
Nonetheless, whenever applicable to the nature of the 
present study, some of the key issues that need to be 
addressed, i.e. preprocessing, preparation and annotation 
of the individual datasets, complied to the guidelines 
suggested by [34]. 



In order to perform data integration, two main issues 
had to be resolved: (i) matching reporters on the two 
microarray platforms, and (ii) normalizing data to 
address platform related differences [35]. 

The first task was rather straightforward, since both 
platforms were obtained from the same manufacturer 
and many reporters were common between the 1.2 k and 
the 4.8 k platforms. Specifically, each reporter-level iden- 
tifier (GenBank accession numbers) was mapped to a 
UniGene identifier (UniGene Cluster ID) [36-40]. The 
mapping was performed through the web-based tool 
SOURCE [41] in the 19th of November 2010, simulta- 
neously for both platforms, in order to avoid inconsisten- 
cies [42]. All mapped reporter-level identifiers had one- 
to-one relationship with the gene-level identifiers that is, 
each reporter was associated with a single UniGene iden- 
tifier and no more than one reporter was mapped to the 
same UniGene identifier. Reporters having insufficient 
information to be mapped to any gene-level identifier 
were omitted. Thus, a fully updated set of unique gene- 
level identifiers was generated for each platform. The 
intersection of these two sets formed a common gene set 
(CGS) consisting of 490 UniGene identifiers common in 
all setups, which was utilized for the subsequent analysis. 

Then, a cross-platform normalization procedure was 
applied in order to address batch effect issues, namely the 
median rank score (MRS) normalization method [43,44]. 
After choosing a microarray set as reference, this simple 
approach replaces the expression values for all the others 
(non-reference) microarray sets by median expression 
values of genes from the reference set. The reference set 
of choice was the 4.8 k platform set, since it outperforms 
in data quality when compared to the 1.2 k one [43]. The 
improved comparability of the data from the two different 
platforms after cross-platform normalization is shown in 
Additional File 1, where it can be seen that the distribu- 
tions of gene expression values derived from the 1.2 k and 
4.8 k platforms are more similar after application of MRS 
in comparison to non-integrated data. 

The following analysis steps were performed on the 
integrated gene expression values of the CGS. 

Identification of differentially expressed genes 

In order to identify potential, differentially expressed 
(DE) genes for each slide, whenever the experimental 
design did not include replicates for the same condition, 
we adjusted a selection process, exploiting the intensity- 
dependent calculation of the standard Z-score [32] of 
each feature within each slide. This approach utilizes a 
sliding window to calculate smoothed local means and 
standard deviations (SDs), which are then used to calcu- 
late the Z-score of the logarithmic ratio values for each 
data point in the normalized MA-plot. At the present 
study, a sliding window of width equal to 20% of the 
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total number of data points was utilized. The selected 
percentage on one hand represents a plausible tradeoff 
between the adequate sample size for statistics calcula- 
tion, and neighbourhood sensitivity regarding the inten- 
sity-dependent action. The DE genes per experiment 
were identified at a confidence level of 95%. 

Cluster analysis 

In order to partition the gene expression profiles 
throughout all five experimental setups and exploit simi- 
larities in the expression profiles, for surmising func- 
tional relevance, possibly through common regulatory 
actions, which orchestrate genomic expression with 
respect to the studied phenotypes in this study, an unsu- 
pervised cluster analysis was utilized. 

The unsupervised clustering was applied to an appro- 
priately selected subset of the CGS, since simulations 
have indicated that keeping irrelevant genes during clus- 
ter analysis, results in reduced accuracy [45] . Specifically, 
each gene of the CGS was ranked by SD across experi- 
ments. The top 100 genes with the highest SD (SD100) 
were selected for downstream analysis. As demonstrated 
in the recent thorough evaluation study for two-channel 
microarray data [46], this widely used method is one of 
the best performing gene selection approaches. 

The method utilized for cluster analysis was the k- 
means clustering [47-49], considered as one of the best 
performing clustering options in particular for microar- 
ray class discovery studies [46]. The k-means algorithm 
applied [50] uses the squared Euclidean as a distance 
measure, which has been frequently found to outper- 
form, regargding ratio-based measurements [51]. Also, 
in order to evade the problem of possible local minima, 
the algorithm is repeated a number of executions with 
different initializations. Specifically, the clustering proce- 
dure is repeated 100 times, each with a new set of initial 
cluster centroid positions (seeds), selected at random, 
and the best execution (the one that minimizes the sum, 
over all clusters, of the within-cluster sums of object-to- 
cluster-centroid distances) is taken as the final result. 

In order to predict the optimal number of clusters for 
interpretation, which constitutes a fundamental problem 
in unsupervised clustering, a cluster validity measure 
was applied. The validity measure of choice was the 
average silhouette width for the entire dataset S com- 
prising N objects [52,53]: 

ieS 

where 5, is the silhouette width of object i, which is 
defined as: 

_ bj - Oj 
1 max(a;, hi) 



where a L is the average distance between object i and 
all the other objects in the cluster, and b t is the mini- 
mum of the average distances between i and objects in 
other clusters. To determine the optimal number of 
clusters [54], the clustering was executed for k varying 
between 2 to 30. For each k, the aforementioned k- 
means algorithm was repeated 1,000 times. The best 
(maximum) values of obtained by each k were 
plotted as the function of k (Additional File 2). As illu- 
strated in this figure, since the plot did not exhibit any 
specific increasing or decreasing trend, we sought for its 
maximum value [54]. Thus, the optimal cluster number 
was the one corresponding to the maximum value of 
the plot and was found to be equal to 7. 

The implementation of all preprocessing steps, along 
with the identification of the DE genes and cluster ana- 
lysis, was performed in the Matlab® 7.9 (R2009b) (The 
MathWorks Inc. Natick, MA, USA) computing 
environment. 

Gene Ontology based analysis 

In order to compare different groups of genes, highlight- 
ing different functionalities among all experimental set- 
ups, interesting genes formed study sets that were 
further subjected to Gene Ontology (GO) based analysis 
to test the nature of the observed resistance mechanism. 
Specifically, for each study set formed, statistical analysis 
of GO term overrepresentation was performed against 
the CGS, which was utilized as a reference set, as pro- 
posed by [55]. The chosen approach was the parent- 
child-union method [56], since it was found to outper- 
form the standard method of overrepresentation analysis 
(ORA) in GO. The standard approach treats each GO 
term independently and hence does not take dependen- 
cies between parent and child terms into account, ignor- 
ing the structure of the GO hierarchy. It was shown that 
this behavior can result in certain types of false positive 
results, with potentially misleading biological interpreta- 
tion [56]. In contrast, the parent-child method measures 
the overrepresentation of a term with respect to the pre- 
sence of its parental terms in the set, and hence resol- 
ving the problem of the standard method which tends 
to falsely detect an overrepresentation of more specific 
terms below terms known to be overrepresented. 

ORA was performed with the publicly available Onto- 
logizer 2.0 tool [57] using GO terms definitions and 
associations between genes and GO downloaded from 
the Gene Ontology consortium [58] on the 26th of 
November 2010. 

TFBMs analysis 

In order to identify the transcription factors driving the 
observed changes in the gene expression, we investi- 
gated the Transcription Factor Binding Motifs (TFBMs) 
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in the Transcription Element Listening System Database 
(TELiS) [59]. The TRANSFAC transcription factor data- 
base was used for the identification of gene transcription 
factor binding sites [60]. 

Chromosome mapping 

Chromosome mapping appears to be a promising 
method for identifying patterns among genes. The main 
idea reported initially by [61] is to map genes in chro- 
mosomal regions and in this way, if correlations do 
exist, they appear through the location of genes on 
chromosomal regions, since consistent expression of 
whole functional entities is related with activation of 
given chromosomal regions. For chromosome mapping 
analyses, we used the Gene Ontology Tree Machine, 
WebGestalt web-tool [62]. 

Pathway analysis 

Pathway analysis was performed in order to find genes 
that participate in known pathways related to the inves- 
tigated mechanisms, such as apoptosis or cell prolifera- 
tion. This approach was used in order to gain further 
insight into the mechanisms underlying the common 
genes identified by our proposed model. The differen- 
tially expressed genes were mapped on different path- 
ways using the Pathway Explorer software [63]. First of 
all, it was investigated the percentage of genes present 
in all known pathways using the databases available 
through the Pathway Explorer software. The KEGG 
Pathways database was used for our analysis [64-66]. 

Hypothesis examination and computational work flow 
Hypothesis statement 

In order to answer the question of inherently disposed or 
acquired resistance and examine our initial working 
hypothesis, significant genes derived from the computa- 
tional workflow were grouped in three sets, for each experi- 
ment i, with i = 1...5, Si, U b and D b where: S t includes genes 
which expression is considered to be unchanged, Zi, 
includes genes up-regulated, and £>, includes genes down- 
regulated. Furthermore, in order to facilitate interpretation, 
gene functions were categorized in two major categories, F p 
and F a , where F p is related to apoptosis evasion and cell 
cycle progression, whereas F a is related to apoptosis induc- 
tion and cell cycle arrest. In other words, F p is related to 
genes that are not desired targets of prednisolone regulation 
and the glucocorticoid receptor (GR) in the system under 
study, while F a is related to genes involved in apoptosis, 
which are the desired ones of GC regulation. 
Modeling approach 

Concerning the 4 h treatment, namely the early treat- 
ment, the following functional categories were defined: 

i) If genes simultaneously remain unaffected at the low 
and high prednisolone doses (intersection of sets Si and 



S 3 (Figure la), which points out dose-independent 
genes), and GO functions are related to F p , then those 
functions seem not to confer to the manifestation of 
resistive behavior by the cells, thus suggesting that these 
cells are originally sensitive and gradually develop resis- 
tant phenotypic characteristics. On the contrary, if gene 
functions are related to F a then cells tend to remain 
unresponsive to GC treatment, thus implying an intrin- 
sic mechanism of resistance. 

ii) If genes are simultaneously up-regulated at both 
concentrations (intersection of sets U\ and U 3 (Figure 
lb), which points out probable early time-dependent 
genes or a generalized reaction mechanism to GCs which 
is dose-independent), and GO functions are related to F p , 
then cells have an intrinsic mechanism of resistance. On 
the contrary, if gene functions are related to F a then cells 
are originally sensitive and gradually develop, or homeos- 
tatically adjust to resistant phenotypic characteristics. 

iii) If genes are simultaneously down-regulated at the low 
and high prednisolone doses (intersection of sets D\ and D 3 
(Figure lc), which points out dose-independent genes, that 
manifest the same expression profile under these two 
extreme GC concentrations, viz, down- regulated in the cell 
system in its natural condition), and GO functions are 
related to F p , then cells are originally sensitive and gradually 
develop resistant phenotypic characteristics. On the con- 
trary, if gene functions are related to F a this predicates 
upon an intrinsic mechanism of resistance by the cells. 

iv) If genes are simultaneously up-regulated at the low 
prednisolone concentration and down-regulated at the 
high prednisolone dose (intersection of sets U\ and D 3 
(Figure Id), which points out the dose-dependent genes) 
and GO functions are related to F p , then cells tend to 
present a dose dependent sensitization to GC exposure, 
thus postponing resistance manifestation at a later point, 
something which implies acquired resistance mechan- 
isms. On the contrary, if gene functions are related to F a 
then cells have an intrinsic mechanism of resistance. 

v) Finally, if genes are simultaneously down- regulated 
at the low prednisolone concentration and up-regulated 
at the high prednisolone dose (intersection of sets Di and 
U 3 (Figure le), which points out the dose-dependent 
genes) and GO functions are related to F p , then cells 
show a dose-dependent intrinsic mechanism of resis- 
tance. On the contrary, if gene functions are related to F a 
then cells present dose-dependent sensitization, postpon- 
ing resistance manifestation at a later point something 
which implies acquired resistance mechanisms. 

Further on, extrapolating the aforementioned 
approach in the analysis of the 72 h experiments, 
namely the late treatment, the following functional cate- 
gories were defined: 

i) If genes remain unchanged (intersection of sets S 4 
and S 5 (Figure If) or are up-regulated (intersection of 
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Figure 1 Venn diagrams illustrating the intersections of the formulations posed for resistance Schematic representation of the set 
intersections (a) - (j) formed based on the definitions described in Methods, in order to answer the question of inherently disposed or acquired 
resistance. Sets 5,-, (J, and D, include genes constantly expressed, up-regulated, and down-regulated in experiment /', respectively. 



sets U 4 and U 5 (Figure lg) at both treatments at 72 h, 
and GO functions are related to F p , then cells manifest 
a time-dependent acquired mechanism of resistance. If 
the opposite is true, then the nature of the resistance 
mechanism cannot be defined. 



ii) If genes are down-regulated (intersection of sets £> 4 
and D 5 (Figure lh) and GO functions are related to F a , 
then the mechanism of resistance is probably inherent. 
Again, if the opposite is true, then the nature of the 
resistance mechanism cannot be defined. 
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iii) If genes are up-regulated at the 22 uM vs. 700 uM 
experiment and down-regulated at the control vs. 700 
uM experiment (intersection of sets U 4 and D 5 (Figure 
li) or are down-regulated at the 22 uM vs. 700 uM 
experiment and up-regulated at the control vs. 700 uM 
experiment (intersection of sets D 4 and ZJ 5 (Figure lj), 
and GO functions are related to F p , then cells manifest 
an acquired resistance mechanism which is partially 
remitting to a dose-dependent re-sensitization by large 
GC doses. On the other hand, if gene functions are 
related to F a , then the mechanism of resistance cannot 
be defined. 

All aforementioned combinations are summarized in 
the simplified flowcharts of Figure 2. 

Results 

Gene expression profiling concerning two distinct time 
points (4 h and 72 h) followed mRNA isolation from 
the CCRF-CEM cell line, cultured in three different pre- 
dnisolone concentrations (10 nM, 22 uM, and 700 uM). 
The mRNA collected was used for hybridization on two 
cDNA microarray platforms (1.2 k and 4.8 k). Specifi- 
cally, the experimental design consisted of the five fol- 
lowing pairs: control vs. 10 nM prednisolone at 4 h 
(designated as '1'), 10 nM vs. 700 uM prednisolone at 4 
h (designated as '2'), control vs. 700 uM prednisolone at 
4 h (designated as '3'), 22 uM vs. 700 uM prednisolone 
at 72 h (designated as '4'), and control vs. 700 uM pre- 
dnisolone at 72 h (designated as '5'). 

A common for both platforms data preprocessing 
stage was applied to the raw expression datasets, includ- 
ing proper background correction for each replicate 
slide, within-slide normalization and filtering. Concern- 
ing data integration, a two-step procedure was applied, 
firstly by matching reporters on the two, microarray 
platforms and then applying a cross-platform, normali- 
zation approach. In this way the preprocessed expres- 
sion values of the common gene-level identifiers, 
existing in both platforms, were selected to form a uni- 
fied dataset on which standard statistical procedures 
could be applied. Each subsequent analysis step was per- 
formed on the integrated gene expression values of this 
unified dataset (Additional file 3), which includes a total 
number of 490 genes, common in all experimental set- 
ups (see Methods for details). 

Gene Ontology analysis based on statistically significant 
genes 

In order to identify differentially expressed (DE) genes 
for each experiment, we utilized the intensity-dependent 
calculation of the standard Z-score [32] of each com- 
mon feature within each slide, at a confidence level of 
95%. Then, for each experiment i, with i = 1...5, all 
genes were further classified into the three sets: S„ £/„ 



and D b where S t includes genes which expression is 
considered to be unchanged, Ui includes genes up-regu- 
lated, and D t includes genes down-regulated (Additional 
file 4). 

Set intersections 

In order to address the question of the prevalence of 
either inherently disposed or acquired mechanisms of 
resistance, the following intersections were defined, in 
line with the formulations described in the Methods sec- 
tion (Figure 1): 5j n S 3 , U 1 n U 3 , D 1 n D 3 , U x n D 3 , D 1 
n U 3 , 5 4 n ,S 5 , U 4 n U 5 , D 4 n D 5 , U 4 n U 5 , and £> 4 n U s 
(Additional file 4). Then, these intersections were sub- 
jected to statistical analysis of Gene Ontology (GO) term 
overrepresentation, in order to test the nature of the 
observed resistance mechanism (Additional file 5). The 
generated from the GO-based analysis gene functions 
were also further categorized in two major categories: F p 
and F a , where F p is related to apoptosis evasion and cell 
cycle progression, whereas F a is related to apoptosis 
induction and cell cycle arrest (Figure 2). 

The results of this examination demonstrated specifi- 
cally that: 

i) Genes belonging to intersection Si n S 3 (Figure la) 
do not incorporate any functions related to apoptosis 
induction (F a ) or cell proliferation (F p ), thus the nature 
of the resistance mechanism cannot be defined. 

ii) No common genes were observed for sets U\ and 
U 3 (Figure lb), implying that no genes from the present 
dataset are dose-independent upon prednisolone 
treatment. 

iii) Intersection D l n D 3 (Figure lc) showed only one 
common gene, the DAPK1. However, DAPK1 is known 
to be a positive mediator of gamma interferon induced 
cell death [67], thus down-regulation of DAPK1 implies a 
possible linkage with inherent resistance mechanisms. 
Moreover, it appears that DAPKl's down-regulation is 
dose-independent and the glucocorticoid receptor (GR), 
when stimulated in the system under study, deactivates it. 

iv) Intersections U\ n D 3 (Figure Id) and D l n U 3 
(Figure le) do not yield common genes. Since there are 
no genes that are regulated dose-dependently and in 
opposing manners, it may be concluded that, at least as 
far as the present dataset is concerned, prednisolone 
regulates different sets of genes in a dose-dependent 
manner. 

Regarding the 72h time point: 

i) Genes belonging to intersection S 4 n S 5 (Figure If) 
comprised functions related to induction of apoptosis 
and cell cycle arrest {F a ), thus they do not confer any 
knowledge about the nature of the resistance 
mechanism. 

ii) Intersections U 4 D Us (Figure lg) and D 4 n D 5 
(Figure lh) did not yield common genes, whereas inter- 
section U 4 n D s (Figure li) captured one gene, namely 



Sifakis et al. Journal of Clinical Bioinformatics 201 1, 1:36 
http://www.jclinbioinformatics.eom/content/1/1/36 



Page 9 of 21 




Figure 2 Flowcharts of the hypotheses posed for resistance. Simplified flowcharts (a) - (d) of the hypotheses posed for inherent or acquired 
resistance in ALL cells. Subsets 5, n S 3 , D, n D 3 and U } n D 3 follow flowchart (a), subsets U ] n U 3 and D, n U 3 follow flowchart (b), subsets S 4 n 
S5, L/ 4 n 1/5, L/4 n D 5 and D 4 n (7 5 follow flowchart (c), and subset D 4 n D 5 follows flowchart (d). Sets S,> L/„ and D, include genes constantly 
expressed, up-regulated, and down-regulated in case /', respectively, where subscript / is varying from 1 to 5, and refers to each experimental 
setup. Each subset is tested for its GO enrichments, while gene functions are categorized as F a , F p , where F a is related to apoptosis induction 
and cell cycle arrest and F p is related to apoptosis evasion and cell cycle proliferation. 



Sifakis et al. Journal of Clinical Bioinformatics 201 1 , 1:36 
http://www.jclinbioinformatics.eom/content/1/1/36 



Page 10 of 21 



the gene KIT. This gene is reported to play a role in a 
variety of human tumors, including acute myelogenous 
leukemia, in its mutated forms [68]. From the present 
study it appears that it is also active in the regulation of 
GC action. 

iii) Finally, one gene is derived from the intersection 
of D 4 n U s (Figure lj), namely MADD (also known as 
IG20), which represents a very interesting gene since it 
has an established role in the mediation of the death 
signal from TNF-alpha downwards to the apoptosis acti- 
vators [69]. In general, it is also known to be expressed 
at higher levels in neoplastic cells than in normal cells 
[70], which is also the case for our system. However, it 
is unclear if MADD over-expression in tumor cells 
incurs increased apoptosis or not. The case of this speci- 
fic gene becomes more interesting as it is known that at 
least six different splice variants with equally different 
functions are expressed, each one performing different 
functions, both apoptotic and anti-apoptotic [70]. 
Individual sets 

From the aforementioned it seems that dose-dependent 
prednisolone administration induces varying expression 
of different sets of genes. In this sense, a plausible pre- 
sumption is that genes differentially expressed by pre- 
dnisolone should to a large extent be implicated in the 
resistance mechanisms, responding to prednisolone 
treatment. Thus, apart from the aforementioned inter- 
sections, the sets U\, U 3 , D 1 and D 4 for the 4 h time 
point, and the sets Z/ 4 , U 5 , D 4 and D 5 for the 72 h time 
point were also examined (Additional file 5). More 
specifically: 

i) Gene set Hi, concerning response to low dose GC 
administration was related to apoptotic and cell cycle 
regulation mechanisms. It is worth noting that this set 
(upregulated by 10 nM prednisolone) includes gene 
BCL2L2, which has been reported to be an anti-apopto- 
tic gene [71]. 

ii) Gene set Z/ 3 , concerning genes overexpressing as a 
result of high GC presence, did not highlighted func- 
tions relevant to apoptosis induction or positive regula- 
tion of cell death (F a ). In this sense, it provides evidence 
for an inherent mechanism of resistance, since it seems 
that the high dose does not comply with the expected 
apoptotic actions. Interestingly, the same concentration 
(700 uM) seemed to stir metabolic as well developmen- 
tal mechanisms. 

iii) Both sets D x and D 3 did not incorporate functions 
related to apoptosis induction and positive cell death 
regulation (F a ). In other words, the genes down-regu- 
lated both in low or high doses of prednisolone adminis- 
tration are not related to apoptosis induction. 
Nevertheless in set Di, gene DAP, which encodes for a 
protein that is positive regulator of programmed cell 
death [72], is related in addition to the process of 



autophagy. Moreover, this gene is a member of the 
mTOR pathway, which is at the same time, a negative 
regulator of autophagy [73]. Again, it appears that 10 
nM prednisolone stimulates cell proliferation in addition 
to an anti-apoptotic effect. It cannot be precluded that 
prednisolone could potentially enhance alternative nutri- 
tion mechanisms as an alternative way to attain survival. 
It appears that GC presence is impacting the way the 
cell will metabolize. 
Regarding the 72h time point: 

i) It appears that prednisolone exerts a diversified 
intricate mechanism as it up-regulates (ii 4 ) genes, such 
as BCL2, which is an anti-apoptotic gene. Moreover, 
functions were revealed that had to do with positive reg- 
ulation of locomotion, indicating the activation of a pos- 
sible mechanism of evasion from a hostile environment. 

ii) Also, prednisolone up-regulated genes related to 
apoptosis evasion (F p ) at the 700 uM dose as compared 
to untreated cells (U 5 ), among these, gene BIRCS. The 
BIRC family of genes belongs to the larger IAP family of 
genes (inhibitors of apoptosis genes). BIRC5 especially is 
considered to be an apoptosis evasion factor with signif- 
icant presence in a variety of tumors [74]. The fact that 
prednisolone up-regulates such a gene, supports the 
case of an acquired mechanism of action. 

iii) In addition, prednisolone down-regulated genes 
{D 4 and D 5 ) are related to stimulation of cell cycle pro- 
gression {F p ), such as BCL2L2 and KIT, respectively. 
Also, genes related to positive induction of cell death, 
such as IKBKE and PPP3CC [40] are down-regulated. 
However, it is worth mentioning that there is also a 
shift of molecular function towards NF-ftB-related 
effects. Especially, gene IKBKE, which is a non-canonical 
Ik-B kinase, a known controller of NF-kB [75], appears 
to be down regulated by prednisolone at 72 h. This is 
the same gene that is involved in induction of cell 
death, i.e. IKBKE, which is a non-canonical Ik-B kinase, 
a known controller of NF-kB [75]. 

A very interesting case is presented in these clusters 
with the MCL1 gene. As mentioned before, this gene is 
instrumentally linked to cell survival and the fact that 
the low dose activates it, is a strong indication for inher- 
ent resistance in the present system. This indication 
becomes stronger by the finding that this gene is simi- 
larly expressed in untreated cells and cells treated with 
various concentrations of prednisolone. As a matter of 
fact it seems that even variations in concentrations such 
as the 22uM and 700uM do not affect its expression, at 
least at late time points, or its expression has been stabi- 
lized. In cluster 7, gene functions appeared to be related 
to cell cycle regulation. This group of clusters seem to 
outline the effects of the high prednisolone dose, as this 
dose activates genes related to cell cycle progression. 
For example, in these clusters appear genes such as 
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BMP5, FGF7 and MCL1. Those genes, are promoters of 
anti-apoptosis and proliferation. Over-expression of 
such genes, at later time intervals, points out the exis- 
tence of an inherent mechanism of reaction to the glu- 
cocorticoid. Yet, the fact that their expression at 4h 
remains stable, implies that for this gene set the decision 
for gene activation or not is taken later on during GC 
treatment. The fact however, is that cells react with 
anti-apoptotic signals already at 4h and our hypothesis 
is that the 72h activation of anti-apoptotic genes has to 
do with a re-enforcement of the resistant phenotype 
that the cells already possess. 

Taken together the results of all the studied subsets, 
summarized in Table 1, showed that the cell system 
under study exerts an intricate response pattern of 
inherent and acquired mechanisms of GC resistance, 
which seems however to favor inherent resistance 
mechanisms. 

We have seen that at early time points, there are no 
genes, which are regulated dose-dependently or in 
opposing manners. Therefore, we concluded that pre- 
dnisolone regulates different sets of genes in a dose- 
dependent manner, at least as far as the present data set 
is concerned. This helped us to conclude that genes dif- 
ferentially expressed by prednisolone should play a role 
in the resistance mechanism to prednisolone. 

Gene Ontology analysis based on clustering 

K-means clustering is an effective way of classifying 
gene expression profiles, based on their similarities as 
well as on probable common regulatory mechanisms. In 
this way, it can be considered as an effective way to sug- 
gest candidate targets of putative, common regulatory 
mechanisms. Thus, as described in Methods, the top 
100 genes with the highest SD (SD100) from the CGS 
was subjected to k-means clustering using squared 
Euclidean as a distance measure, in order to identify 
groups of genes presenting similar expression profiles 
during prednisolone treatment early and late response 
and possibly co-regulated. In Figure 3, the k-means clus- 
tering is presented, illustrating the seven clusters derived 
(Additional file 6). GO-based enrichment analysis, that 
followed, linked co-expressed genes in each cluster to 
several functional categories (Additional file 5). In parti- 
cular, we noted that: 

i) Cluster 1 (Ci) contains time-dependent genes, which 
expression decreases at 700 uM prednisolone and the 72 
h time point (resembling the behaviour of the D 5 set), 
while remaining unchanged at all other conditions 
(including genes from the Si n S 2 set). However, genes 
under this cluster did not yield any significant GO terms. 

ii) Cluster 2 (C2) mainly comprises genes showing 
early down-regulation at 700 uM prednisolone (includ- 
ing genes from the D 3 set), while remaining unchanged 



Table 1 Summary of the GO functions revealed for each 
gene set 

Gene Group 3 No. Gene b Gene Function 1 Resistance* 1 p-value e 



* ~a f"m 'i 



4 h S, n 

S3 


423 




no 


no 


no 


yes 


n/d 


< 0.05 


Uj n 


0 




n/ 


n/ 


n/ 


n/ 


n/a 


n/a 


U3 






a 


a 


a 


a 






D, n 

D 3 


1 


DAPK1 


no 


yes 


no 


no 


inherent 


< 0.05 


U, n 


0 




n/ 


n/ 


n/ 


IV' 


n/a 


n/a 


D 3 






a 


a 


a 


a 






D, n 


0 




n/ 


n/ 


n/ 


n/ 


n/a 


n/a 


U 3 






a 


a 


a 


a 






u, 


12 


BCL2L2 


yes 


no 


no 


no 


inherent 


< 0.05 


u 3 


11 




yes 


no 


yes 


no 


inherent 


< 0.05 


D, 


23 


DAP 


no 


yes 


yes 


no 


inherent 


< 0.05 






CASP1 


no 


yes 


no 


no 


inherent 


< 0.05 






CDC25 


yes 


no 


no 


no 


acquired 


< 0.05 


D 3 


22 


1124 


no 


yes 


no 


no 


inherent 


< 0.05 


72 h S 4 n 

5 5 


433 




no 


yes 


yes 


no 


n/d 


< 0.05 


U 4 n 


0 




n/ 


n/ 


n/ 


n/ 


n/a 


n/a 


Us 






a 


a 


a 


a 






D 4 n 

D 5 


0 




n/ 
a 


n/ 
a 


n/ 
a 


IV' 

a 


n/a 


n/a 


U 4 n 

D 5 


1 


KIT 


yes 


no 


no 


no 


acquired 


< 0.05 


D 4 n 

u i 


1 


MADD 


n/ 

A 
U 


n/ 

(J 


n/ 

(J 


n/ 

(J 


n/cl 


< 0.05 


U 4 


9 


BCL2 


yes 


no 


no 


no 


acquired 


< 0.05 


Us 


12 


BIRC5 


yes 


no 


no 


no 


acquired 


< 0.05 


D 4 


19 


BCL2L2 


yes 


no 


no 


no 


acquired 


< 0.05 


Ds 


18 


IKBKE 


no 


yes 


no 


yes 


~ inherent 


< 0.05 






PPP3CC 


no 


yes 


no 


no 


~ inherent 


< 0.05 


clusters c l 


A 




n/ 
a 


n/ 
a 


n/ 
a 


n/ 
a 


n/a 


n/a 


c 2 


11 




n/ 
d 


n/ 

cl 


n/ 
d 


n/ 
d 


n/a 


n/a 


c 3 


14 




no 


yes 


no 


no 


inherent 


< 0.05 


c 4 


25 




n/ 
a 


n/ 
a 


n/ 
a 


n/ 
a 


n/a 


n/a 


c 5 


17 




yes 


yes 


no 


no 


n/cl 


< 0.05 


c 6 


10 




no 


yes 


no 


no 


inherent 


< 0.05 


c 7 


19 




yes 


no 


yes 


no 


acquired 


< 0.05 



a Sj, Uj, and D, correspond to groups which include genes constantly 
expressed, up-regulated, and down-regulated in experiment /", respectively, 
and Cj to a cluster of genes showing similar expression profile. Subscript / is 
varying from 1 to 5 and refers to each experimental setup, b Selected genes 
from each gene set, c F p = is related to apoptosis evasion and cell cycle 
progression, F a = is related to apoptosis induction and cell cycle arrest, F m = is 
related to metabolic processes, F-, = is related to the regulation of IkB/NFk.B, 
n/a = not available, and n/d = not defined, d ~ inherent = probably inherent, 
e Each p-value corresponds either to the one-sided Fisher exact test of the 
GO-based analysis of the respective gene set, or, in case of a single gene, to 
the intensity-dependent Z-score. 
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1 2 3 4 5 

Figure 3 Expression profile clustering. K-means clustering of the gene expression profiles illustrating the seven clusters as described in 
Methods. Each subplot corresponds to one cluster of genes showing similar expression profile. The horizontal axes correspond to the five 
prednisolone experiments of interest, where T = control vs. 10 nM prednisolone at 4 h, '2' = 10 nM vs. 700 uM prednisolone at 4 h, '3' = 
control vs. 700 uM prednisolone at 4 h, '4' = 22 uM vs. 700 uM prednisolone at 72 h, and '5' = control vs. 700 uM prednisolone at 72 h. The 
vertical axes depict the corresponding logarithmic ratios, as derived from microarray data analysis. 



at all other conditions (probably outlining genes of the 
S 4 n S s set). GO-based analysis revealed genes involved 
in developmental processes, implying a possible exis- 
tence of stem cell related functions in the cell system 
under study. 

hi) Cluster 3 (C 3 ) depicts genes showing early down- 
regulation at 10 nM prednisolone (including genes from 
the Di set), while remaining unaffected by the 700 uM 
prednisolone treatment. This cluster comprised genes 
related to positive regulation of cell death. The fact that 
these genes are suppressed under low prednisolone con- 
centration at the early time point is in agreement with 
the anti-apoptotic and survival effect reported previously 
by [76], which implies an intrinsic mechanism of 
resistance. 

iv) Cluster 4 (C 4 ), similarly to (Q), presents genes 
which expression at 700 uM prednisolone decreases at 
the 72 h time point (resembling the behaviour of the D s 
subset), while remaining unchanged at all other condi- 
tions (including genes from the Si n S 3 subset). How- 
ever, genes under this cluster did not present any 
interesting GO terms. 



v) Cluster 5 (C 5 ) contains genes remaining unchanged 
at all conditions (including genes from the Si n S 3 inter- 
section) except for experimental setup '4', where they 
are down-regulated (resembling the behaviour of the £> 4 
set). GO-based analysis of cluster 5 revealed genes that 
are involved in the opposing processes of cell cycle pro- 
gression, as well as cell cycle arrest, and cell death 
functions. 

vi) Cluster 6 (Cg) depicts genes that are also 
unchanged at all conditions (including genes from the 
5 4 n S 5 intersection) except for experimental setup '2', 
where they seem to be down-regulated (resembling the 
behaviour of the D 2 set). Functions represented by genes 
in cluster 6 include cell death functions so their down 
regulation at the specific setup supports the case of 
inherent resistance mechanisms. 

vii) Finally, cluster 7 (C 7 ) groups together time-depen- 
dent genes that remain unaffected at all experimental 
setups (including genes from the Si n S3 intersection), 
but at 700 uM after 72 h of treatment (resembling the 
behaviour of the U 5 set), where they are positively regu- 
lated. It appears that those genes participate in cell cycle 
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regulation. This cluster seems to outline the effects of 
the high prednisolone dose, as this dose activates genes 
related to cell cycle progression. Special attention was 
drawn to the MCL1 gene. This gene is a member of the 
BCL2 family and it produces an anti-apoptotic protein 
responsible for cell survival. 

The cluster analysis results, summarized in Table 1, 
also confirm that, although the gene expression profiles 
of the cell system under study present a mixed response 
of inherent and acquired mechanisms of GC resistance, 
however there is a tendency towards inherent resistance 
mechanisms. 

TFBM analysis 

A next step on our analysis was the identification of 
common transcription factors that would regulate genes 
in a similar way. Cluster 1 (Ci) manifested a common 
transcription factor namely OCT1 (p = 5 • 10" ). Cluster 
2 (C 2 ) similarly manifested a common transcription fac- 
tor, MZF1 [p = 8 • 10" 4 ). Cluster 5 (C 5 ) were predicted 
to be commonly regulated by the CCAAT (NFYA) tran- 
scription factor (p = 5 • 10~ 7 ). Cluster 6 (C 6 ) presented 
an interesting case, as the common most prevalent tran- 
scription factor was JUN. Finally, Cluster 7 (C 7 ) also 
manifested different transcription factors. In particular, 
the most prevalent transcription factor was OCT1 (p = 
1 • 10' 10 ). 

Chromosome mapping 

In order to find further patterns of expression in the 
system under study, we have performed a chromosome 
mapping analysis, where the mean expression of genes 
with respect to chromosomes has been studied. A gen- 
eral remark from this analysis is that expression seemed 
to be equivalently distributed; at least as far as positive 
regulation is concerned, while chromosomes 14 and 22 
appeared to have genes with most deactivation as com- 
pared to untreated genes (Figure 4a). 

Further on, we have separated gene expression based 
on the time points i.e. we have calculated the mean 
expressions for the 4 h and 72 h treatments (Figure 4b 
and 4c, respectively). We have tried to identify whether 
the mean expression values, either positive or negative, 
correlate to the number of genes represented on each 
chromosome. In Figure 4d, the number of genes per 
chromosome is presented. Table 2 presents the calcu- 
lated correlation values of mean expression values for 
experimental time points, while Figure 5 graphically 
illustrates the Pearson's correlation values comparing 
gene numbers to mean gene expression. Interestingly, 
gene number is negatively correlated to negative mean 
expression, which means that the more genes are repre- 
sented on chromosomes the more these genes are 
down-regulated as compared to untreated cells or, in 



order to be more concise, the more these genes are 
inactivated as compared to control. 

Pathway analysis 

In order to gain further insight into the possible 
mechanisms that probe for the preponderance of inher- 
ence or acquired resistance mechanisms, we have 
attempted to perform a pathway analysis with the 
defined genes from our dataset. We have mapped the 
common dataset on known pathways and in particular 
on the KEGG Pathways database. In particular, we have 
outlined two pathways: the MAPK pathway, as it incor- 
porated the larger number of involved genes to be 
mapped on it (Figure 6), and the apoptosis pathway, 
since we were interested in examining the possible role 
of genes related to apoptotic mechanisms (Figure 7). 

It is known that the MAPK signaling pathway is 
involved in cell cycle progression and apoptosis. We 
have outlined several cross-talking molecules in this 
pathway that were revealed also by our dataset, such as 
the NF-kB, JNK, p38 and ERK5. NF-kB remains stable 
across the experimental setups indicating a steady role 
in apoptosis. Yet, as it is regulated by the GR then its 
levels should be lowering as concentration increases. 
This hints towards a more specific, more pronounced 
involvement of the NF-kB factor in the observed resis- 
tance. At the same time, the JNK kinase, which is a 
mediator of apoptosis, appears to remain unaffected by 
the prednisolone concentrations indicating that the 
MAPK pathway involved in apoptosis does not function 
in the way it should have, hence supporting an inherent 
mechanism of resistance. Finally, an interesting molecule 
is the ERK5, which appears to be down-regulated at the 
low doses and unaffected thereafter. This molecule is 
involved in cell proliferation and its down-regulation is 
one of the very few implications that prednisolone 
behaves as expected, hence indicating an acquired 
mechanism of resistance. 

Moving to the apoptosis pathway, several interesting 
genes that play a role in the observed apoptosis were 
identified. In particular, based on the signaling on Figure 
7, we can discriminate between the TNF/FADD/CASP/ 
BIRC system (top of the figure), where it appears that 
molecules responsible for degradation are down-regu- 
lated under prednisolone treatment. In particular, BIRC5 
was found to be activated in the 72 h treatment, while 
remained unchanged in the previous time point. In 
addition, CASP1 was found to be down-regulated at 10 
nM treatment and 4 h, indicating an inherent mechan- 
ism. This part of the pathway seems to function in such 
a way as cells already possess glucocorticoid resistance. 
The second part of the pathway consisting from NGBF/ 
NTRK1/AKT1/NFkB/BIRC/BCL2 (family) shows the 
following. AKT1 appeared to remain unchanged at all 
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Figure 4 Chromosome-related gene expression. The mean values of all genes at all time points (a) is presented separated in positive and 
negative fold expression. Chromosomes 14 and 22 showed to have the most deactivated genes as compared to untreated cells. Similarly, mean 
gene expression has been separated to the 4 h (b) and 72 h treatments (c). Further on we have searched whether mean gene expression is 
associated with the represented genes on each chromosome as the number of genes is presented on (d). 



concentrations and time points. NF-kB1 is down-regu- 
lated at the late time points, indicating a late response 
to the glucocorticoid, while it is unchanged at the early 
time points. The BCL family is a super-family of pro- 
teins that when expressed promotes survival. In the 



present dataset, the BCL2 gene was found to remain 
unchanged, while a member of the super-family, 
BCL2L2, was up-regulated. BCL2L2 is one of the genes 
that participate in survival in the pathway, confirming 
the fact that we have an early anti-apoptosis response. 
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Table 2 Pearson's correlation analysis of the mean expression values of genes with respect to their chromosome 
allocation 

Pearson's correlation value 3 

All 4 h 72 h 







Genes no. 


m M > 0 


m M < 0 


m M > 0 


m M < 0 


m M > 0 


m M < 0 




Genes no 


1 .0000 


0.2813 


-0.0368 


0.0719 


0.2309 


0.3883 


-0.4294 


All 


m M > 0 


0.2813 


1 .0000 


0.0383 


0.4537 


0.1325 


0.7978 


-0.2066 




m M < 0 


-0.0368 


0.0383 


1 .0000 


0.2503 


0.7349 


-0.0249 


0.5293 


4 h 


m M > 0 


0.0719 


0.4537 


0.2503 


1 .0000 


0.2789 


0.0053 


-0.0059 




m M < 0 


0.2309 


0.1325 


0.7349 


0.2789 


1 .0000 


0.0385 


-0.0599 


72 h 


m M > 0 


0.3883 


0.7978 


-0.0249 


0.0053 


0.0385 


1 .0000 


-0.2422 




m M < 0 


-0.4294 


-0.2066 


0.5293 


-0.0059 


-0.0599 


-0.2422 


1.0000 










p-value a 


















All 




4 h 




72 h 






Genes no 


m M > 0 


m M < 0 


m M > 0 


m M < 0 


m M > 0 


m M < 0 




Genes no 


0.0000 


0.1936 


0.8677 


0.7444 


0.2891 


0.0671 


0.0409 


All 


m M > 0 


0.1936 


0.0000 


0.8624 


0.0297 


0.5467 


0.0000 


0.3441 




m M < 0 


0.8677 


0.8624 


0.0000 


0.2493 


0.0001 


0.9101 


0.0094 


4 h 


m M > 0 


0.7444 


0.0297 


0.2493 


0.0000 


0.1975 


0.9808 


0.9786 




m M < 0 


0.2891 


0.5467 


0.0001 


0.1975 


0.0000 


0.8616 


0.7859 


72 h 


m M > 0 


0.0671 


0.0000 


0.9101 


0.9808 


0.8616 


0.0000 


0.2655 




m M < 0 


0.0409 


0.3441 


0.0094 


0.9786 


0.7859 


0.2655 


0.0000 



a m M = mean gene expression. 



The role of late NFKB1 inactivation and early BCL2L2 
activation probably favors an acquired mechanism of 
resistance. Finally, an apoptotic signaling avenue is via 
the regulation of Ca +2 where the caspases participating 
in this pathway are inactivated in the present system. 
The PPP3CC protein that regulates phosphorylation of 
down-stream targets has a relation to apoptosis over the 




Figure 5 Pearson's correlation values comparing gene 
numbers to mean gene expression Graphical illustration of 
Pearson's correlation values between gene numbers represented on 
each chromosome and the mean expression values for 
experimental setups. This analysis showed no significant results 
except for the case where negative gene expression is correlated to 
gene number. That means that the more genes are represented on 
a chromosome the less these are expressed. 



activation of CASP proteins. Both are inactivated in the 
present system with PPP3CC being inactivated at 72 h 
and CASP1 at 4 h. This presents a mixed mechanism of 
apoptosis evasion. 

Discussion 

In the present work, we have set up and propose a 
rational computational analysis framework, in order to 
aid the elucidation of the molecular mechanisms of glu- 
cocorticoid resistance and specifically whether these are 
to a larger extent inherent or acquired. To the best of 
our knowledge, there are no reports trying to analyze 
systematically the underlying, inherent or acquired 
molecular mechanisms of GC resistance. The issue 
whether leukemic cells possess the inherent genetically 
imprinted resistance mechanisms or it is rather a post 
effect of evolutionary adoption of originally sensitive 
cells upon GC treatment is still a controversial one, 
with potentially significant interest for design of novel 
therapeutic approaches in cancer treatment. A similar 
work [77] reported recently that in ex vivo samples, leu- 
kemic cells exhibit an at least in part, intrinsic mechan- 
ism of reaction. In another work [78], it has also been 
mentioned that glucocorticoids can induce intrinsic 
mechanisms of GC-resistance such as the BCL2 relais. 

The initial working hypothesis was that the in vitro 
system of the present study presents an inherent 
mechanism of reaction to glucocorticoids, resulting in 
resistance to apoptosis. Our results showed that cells 
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Figure 6 Mapped genes of the common dataset on the KEGG MAPK pathway. 



may attain a mixed mechanism of response to glucocor- 
ticoids, however, there is clear evidence predicating 
towards an intrinsic mechanism of resistance. 

Specifically, several interesting genes, whose expres- 
sion is regulated by the GR, were identified. One of 
these genes was MCL1. This gene belongs to the BCL 
super-family, which is responsible for cell survival. It has 
been reported that down-regulation of MCL1 sensitizes 
T-cell leukemia cells to treatment with glucocorticoids 
[79]. In our system this gene appeared to be stable 
across all concentrations while, interestingly, it was up- 
regulated by the low prednisolone dose, confirming the 
anti-apoptotic effect observed previously by the low 
dose. At the same time another member of the BCL 
family, BCL2L2, also responsible for cell survival was 
up-regulated by the low dose of prednisolone. There are 
no reports connecting BCL2L2 with leukemia, yet some 
reports refer to its pro-survival role in leukemogenesis 
[80,81]. An interesting gene also revealed by the present 
analysis was DAPK1, a tumor suppressor gene, 



simultaneously down-reglulated at both prednisolone 
concentrations. Methylations of this gene have been 
linked to hereditary character of chronic leukemias 
[82,83] as well as its expression has been linked to child- 
hood leukemia [84,85]. Additionally, when moving to 72 
h of treatment a gene down-regulated by prednisolone 
was KIT, a proto-oncogene homolog to c-KIT. Its inhibi- 
tion has been reported to be involved in leukemia treat- 
ment and glucocorticoid activity [86-88]. Further on, a 
gene that is up-regulated by the high prednisolone dose 
was MADD. This gene is a propagating agent of the 
death signal induced by TNF. It mediates the signal 
from TNF to MAPK pathway thus inducing apoptosis. 
It is reported to be overexpressed in several neoplasms 
[89] and also it is reported that its over-expression is 
linked to tumor survival [69]. 

MADD is a propagating agent of TNF induced death 
signals, and more specifically, it mediates the signal 
from TNF to MAPK pathway, thus inducing apoptosis. 
MADD is phosphorylated by AKT which then inhibits 
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TRAIL-induced apoptosis [69]. Considering the fact that 
MADD's expression pattern is in line with the resistance 
observed in this study, a lingering question is if this 
over-expression is the result of an inherent or acquired 
mechanism. 

In the search for shared transcription factor binding 
motifs an interesting finding was that of JUN. This fac- 
tor has been reported previously to play a role in the 
apoptosis regulation in the same system studied here 
[90]. 

Chromosome mapping and chromosomal-relative 
expression did not correlate with the number of genes 
represented on each chromosome, thus suggesting a 
well-coordinated GR-induced regulation. GCs are 
known to lead to chromatin modification about 2 h fol- 
lowing treatment. A detectable change in the expression 
of genes that are immediately regulated by the GC is 
expected to reach the stage of mRNA steady-state levels 
at approximately 3-4 h after treatment. Then, genes 
regulated by the GC, begin to influence subsequent gene 



expression. cDNA microarray analysis upon 4 h of treat- 
ment has been chosen, in order to reduce the complex- 
ities that arise later due to ensuing feedback 
mechanisms, at the same time focusing on the GR/NF- 
/tB-linked, direct target genes [76]. Moreover, the delay 
in the action of GCs has been reported depending on 
the cell type and lasting from 2 h to 24 h [91]. Specifi- 
cally, in our model a prednisolone action lag lasts up to 
72 h [76], as reported in different studies [92]. Chromo- 
somal-related expression, is closely linked to chromatin 
remodeling and it could be the object of future investi- 
gations, as it appears that the GR affects almost all chro- 
mosomes at all time points as could be easily 
hypothesized from the abundance of GR response ele- 
ments throughout the human genome [10]. This con- 
firms the idea of complete genome regulation by the 
glucocorticoid receptor in cellular systems. 

The final remarks regarding resistance to glucocorti- 
coids, comes from pathway analysis. Based on our 
observations, probably the apoptosis induced by either 
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TNF receptors or NGFB over the AKT pathway and in 
concordance to the MAPK pathway is ruled out, since 
AKT remains unchanged. That means that in the pre- 
sent system, glucocorticoids do not significantly affect 
this pathway. It is also reported that protein-kinase net- 
works control in a major way the observed resistance to 
glucocorticoids [93]. On the other hand, based on the 
pathway model, two alternative ways remain; the mito- 
chondria-induced apoptotic pathway and Ca +2 induced 
pathway. From those two the mitochondria-directed 
pathway appears to entail genes that remain unchanged 
from glucocorticoid treatment. If we also take into 
account the notion that GR translocation to the mito- 
chondria is part of the resistance mechanisms to gluco- 
corticoids [93] then this action could be related to the 
induction of resistance in our system, as through trans- 
location to mitochondria GR ceases to exert its regula- 
tory effects. Alternatively, a likely mechanism for 
induction of resistance to glucocorticoid-induced apop- 
tosis could be the Ca +2 signaling pathway. Several key 
genes in it appear to be differentially expressed, such as 
PPP3CC and CASP1, at the early time points, thus mak- 
ing evasion of apoptosis probable. 

Attention should be also drawn on two categories of 
genes regulated by prednisolone. These are metabolic 
genes and signal-transduction related genes. In both 
time points, high prednisolone concentration regulates 
such genes, thus grounding for cell proliferation 
machinery. In addition, the regulation of NF-ftB-related 
genes implies an inherent mechanism of resistance 
through the established link of NF-kB inflammatory role 
and GC-induced resistance. 

Overall, the results of the present study support our 
initial working hypothesis for a rather inherent GC 
resistance mechanism. However, further exploitation of 
the proposed computational analysis workflow through 
the conduct of a larger genome-wide transcriptomic 
study could promote the extent and level of detail 
regarding our knowledge about the molecular underpin- 
nings which orchestrate the manifestation and dynamics 
of the mechanisms of GC resistance and desensitization. 

Conclusions 

The leukemic cells used in this study are known to be 
resistant to glucocorticoids and in particular to predniso- 
lone. In order to gain more insight to the mechanisms of 
resistance we have developed a rational computational 
analysis framework along with experimental approaches 
in order to answer the question of whether cells exhibit 
this resistance due to an inherent or, an over time, 
acquired mechanism. We have used the biological infor- 
mation derived from our modeling approach to interpret 
the findings observed regarding our initial hypothesis. 
The analysis of the results supports a complex 



mechanism of action for the cells which tends to favor an 
intrinsic mechanism of resistance, although in one case 
prednisolone appeared to exert the expected effects i.e. 
positive regulation of apoptotic genes, thus supporting an 
acquired mechanism. It seems that a crucial determinant 
for the manifestation of either phenotype seems to be the 
dosage as other haphazard environmental factors, which 
partly regulate the cellular circuitry towards either direc- 
tion. The ability to discriminate between acquired or 
intrinsic mechanisms of resistance is of major importance 
both in the mechanics of glucocorticoid signal transduc- 
tion as well as to the clinical praxis since GCs are still 
front line medications for the treatment of malignant dis- 
eases, especially in the case of leukemia. 

Additional material 



Additional file 1: Cross-platform normalization One microarray slide 
per platform (1.2 k and 4.8 k) was selected and a quantile-quantile plot 
(QQ-plot) was produced (a) before and (b) after the application of cross- 
platform normalization. In each QQ-plot, the quantiles of all gene 
expression values of the first slide were plotted against the quantiles of 
all gene expression values of the second slide. In the case where the 
gene expression values of the two slides come from the same 
distribution, the points in the plot should fall near the straight line. 

Additional file 2: Optimal cluster number determination. K -means 
clustering was executed for a number of clusters, varying between 2 to 
30. For each cluster number, the best (maximum) value of all the 
average silhouette widths obtained at 1,000 executions was plotted 
against the cluster number. Since the maximum values of the average 
silhouette width did not exhibit any specific trend, the optimal cluster 
number was determined as the one corresponding to the maximum 
value of the plot, indicated by the arrow. For the computation of the 
silhouettes the squared Euclidean distance was also used. 

Additional file 3: Unified dataset after data integration The list of 
the 490 genes, common in all experimental platforms and replicates 
(CGS) and their corresponding fold change ratios per experiment. Data 
integration was performed after (i) matching the reporters on the two 
microarray platforms through UniGene Cluster IDs, and (ii) applying a 
cross-platform normalization approach. 

Additional file 4: Gene subsets based on the formulations posed for 
resistance. The list of genes per subset as derived after intensity- 
dependent calculation of the standard Z-score, along with the 
intersections based on the formulations described in Methods. 

Additional file 5: GO terms predicted for each gene subset The list 
of the most significant (p < 0.05) GO terms per gene subset, as derived 
from the GO enrichment analysis. 

Additional file 6: Cluster membership of the top 100 genes with 
the highest standard deviation. The list of the SD100 gene set, along 
with the list of genes per cluster showing similar expression profile 
according to cluster analysis. 
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